function [psi,psi_x,psi_y] = psi_h(x,y,Re,caseNum)
% function [psi,psi_x,psi_y] = psi_h(x,y,Re,caseNum)
%
%
%%%%%%%%%%%%%%%%%%%%%%%%
    
    switch caseNum
        case 0    
            % analytic solution in paper of Shao
            % the domain is [0,1]x[0,1]
 %           nu = 1/Re;
            val = x.*y.*(x-1).*(y-1);
            psi = val.*val;
            psi_x = 2*val.*y.*(y-1).*(2*x-1);
            psi_y = 2*val.*x.*(x-1).*(2*y-1);
%           psi = x.*x.*(x-1).*(x-1).*y.*y.*(y-1).*(y-1);
%           psi_x = 2*x.*y.*y.*(x - 1).*(y-1).*(y-1).*(2*x-1);
%           psi_y = 2*x.*x.*y.*(x - 1).*(x-1).*(y-1).*(2*y-1);
        case 1   
            % analytic solution of Kovasznay
            % the domain is [-0.5,1.5]x[0,2]
            lambda = Re*0.5-sqrt(Re*Re*0.25 + 4*pi*pi);
            psi = y - 0.5/pi*exp(lambda*x).*sin(2*pi*y);
            psi_x = x*0;
            psi_y = y*0;
        otherwise
            psi = 0*x;
            psi_x = x*0;
            psi_y = y*0;
    end
end